function diff = J_r_p(m,param)

    % Derivative of J in the replacement region

    diff = ((1-param.omega_1)./rho_fun(1,'R',param) ...
            + param.J_rep_1*m.^(param.gamma_R_1).*param.gamma_R_1./m + param.J_rep_2*m.^(param.gamma_R_2).*param.gamma_R_2./m) ...
            +J_0_p_fun(m,param);

end